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Abstract 

An in-depth analytic study of a model of language dynamics is 
presented: a model which tackles the problem of the coexistence of 
two languages within a closed community of speakers taking into ac- 
count bilingualism and incorporating a parameter to measure the dis- 
tance between languages. After previous numerical simulations, the 
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model yielded that coexistence might lead to survival of both lan- 
guages within monolingual speakers along with a bilingual community 
or to extinction of the weakest tongue depending on different param- 
eters. In this paper, such study is closed with thorough analytical 
calculations to settle the results in a robust way and previous results 
are refined with some modifications. From the present analysis it is 
possible to almost completely assay the number and nature of the equi- 
librium points of the model, which depend on its parameters, as well 
as to build a phase space based on them. Also, we obtain conclusions 
on the way the languages evolve with time. Our rigorous consider- 
ations also suggest ways to further improve the model and facilitate 
the comparison of its consequences with those from other approaches 
or with real data. 

1 Introduction 

At whatever scale that we look, languages reveal themselves as very elabo- 
rated entities consisting of many coupled parts: grammar, vocabulary, etc; 
each of them complex in its own nature as well. They are, moreover, a main 
instrument of interaction in an entangled web of social agents so that the 
state and evolution of tongues cannot, ultimately, be considered as detached 
from other social dynamics. We readily appreciate that we are in front of an 
utter challenge to the human intellect [TJ [2] . Small steps are gradually taken 
towards a further understanding of the many problems posed by languages. 
Leaving aside those contributions from the more classic fields (e.g. philol- 
ogy), linguistic questions were opened to very diverse branches of science 
during the 20th century by drawing inspiration from some pioneer multi- 
disciplinary works [3]. Given the complexity outlined before, any of these 
first transversal approaches are necessarily simplistic or rely largely on com- 
puter simulations, and rigorous and definitive mathematical proofs of the 
results are often missing. 

The kind of questions that were exposed to a more varied community of 
researchers regard the evolution of languages: transformations in their syn- 
taxes, grammars, or vocabularies; aging, rise, and death; the dynamics of 
their number of speakers: spreading of culture, competition or other kind of 
interaction with other tongues; etc. And the fields that take on these issues 
are as diverse as sociology, biology, or physics. The science of complex sys- 
tems should be highlighted because of its very clever usage of existing mathe- 
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matical methods that stem mainly from statistical mechanics [4]. While also 
appraising other important contributions that widen our knowledge about 
the nature of human languages [HI El El IH1 El Ell EEH EE2|, we shall focus on 
a seminal paper by Abrams and Strogatz [13j that prominently triggered re- 
search in its direction. For an exhaustive review on very varied related topics 
with up-to-date bibliography consult [14]; and for a more extended review 
on the impact of statistical physics on social dynamics, including language 
modeling, see [T5] . 

The line of research propelled by [13] addresses the modeling of language 
coexistence as a competitive dynamics to attract speakers. In [13] a minimal 
model was accounted for and in accordance with experimental data a sounded 
result spread: that a two-languages competition for speakers always led to the 
extinction of one of the parties. Further analysis of the model [HI |T7] shows 
that it also allows for stable language coexistence, but the parametric setup 
needed has not been observed in any study with available real data. Following 
the trend, more complicated models were developed that took spatial or social 
structure into account [THl ESI ED] or that explicitly introduced bilingualism 
|2Ql ED 1221 E31 EU E51 EE]- This naturally eased the way to solutions with 
stable coexisting languages. 

As it was advanced before, computer simulations constitute a favorite tool 
in this modern wave of scientific approaches to the study of languages. The 
results are usually convincing more than enough and, besides, these numerical 
studies allow to reach a depth of knowledge that might be impossible if we 
should rely only on very rigorous analytical demonstrations. Despite of this, 
many of the most insightful contributions to the comprehension of human 
communication follow from meticulous and carefully proven mathematical 
constructions, mainly in the study of grammars and largely aided by methods 
from computational sciences [27] . 

In this paper we intend to make a contribution by analytically elucidating 
some existing results in the modeling of language competition. This is nec- 
essarily a rearguard job-since computer simulations have the lead by far-but 
we will see how it is very valuable and necessary. Thanks to the thorough 
reasoning of this paper we gain a deep understanding about the dynamics of 
speakers of coexisting languages. We focus on a model of language competi- 
tion that allows bilingualism introduced in [23] and whose most interesting 
results were numerically derived in ]25] . By analyzing the model we will come 
to a better interpretation of the previous numerical work and we will reach 
some new results that are, now, supported by robust analytical proofs. 



3 



The paper is structured in the following way: In section [2] the hypothe- 
sis used in [23] are formulated and the corresponding equations are derived 
therefrom. In section [3] strict mathematical results are carefully obtained. 
Whenever the analytic tools do not reach to fully solve the problem, numer- 
ical simulations are employed, but its presence-usually at the very end of 
the chain of reasoning-is always warned to the reader to leave any minimally 
uncertain result open to debate. In section H] the more mathematical aspects 
of the reached solutions are left aside and the results are analyzed primar- 
ily from the point of view of language dynamics: what does the analytical 
outcome mean in terms of coexisting languages? 

2 Derivation of the model equations 

Closely following the path pointed out by Abrams and Strogatz [13], we work 
on a model of two competing languages where bilingualism is an option in 
between and where the similarity between the tongues plays and explicit role 
[23] . The model considers a human population whose individuals might talk 
either of two languages X or Y or both of them. Along the text we might 
refer to either of the monolingual communities or the bilingual one as groups, 
standing for groups of speakers. Naming x and y the fraction of monolingual 
speakers of each tongue, and naming h the fraction of bilingual speakers; two 
non- linear coupled differential equations are derived from the following basic 
hypothesis: 

1. Population size remains constant. 

2. The probability that an individual acquires a language different from its 
current one grows with the status of the new language. Therefore a 
status parameter s G [0, 1] of one of the languages is introduced (being 
1 — s the status of the other one). These statuses are constant and a 
property of the system of coexisting languages. 

3. It is possible to define a distance between two languages. This was 
done in |23j introducing a parameter called interlinguistic similarity, 
k e [0, 1]: k = for orthogonal languages, k = 1 for exactly equal 
languages-a measure of how close the languages are to each other. This 
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allowed to get this distance in an easy and straight way, by simply fit- 
ting percentages of speakers of the involved languages along time to a 
set of differential equations. This interlinguistic similarity is constant 
and, again, a property of each pair of languages. This parameter de- 
scribes how difficult it is for a monolingual speaker to learn the other 
language: this task should be easier if the languages are more similar 
to each other. Another view of it is that the probability that an in- 
dividual retains its old language when learning a new one grows with 
this parameter k. Thus, the parameter k is the gate which opens the 
way to the birth of a bilingual group. 

4. The probability that an individual acquires a language different from its 
current one grows with the fraction of speakers of the new language. 
An exponent a is introduced to ponder the importance of the fraction 
of speakers in a group in attracting new speakers with respect to s and 
k. Once more, a should be a constant that characterizes each system 
of two competing languages. Existing field work with similar equations 
[T5] shows that the equivalent parameter in that model varies little 
across pairs of coexisting languages, suggesting that social pressure to 
shift languages might be a constant through cultures. 

Considering different hypothesis or variations on the implementation of 
the current ones might lead to different modeling of the same phenomenon 
[2TI l22l [23]. In the present paper we focus on this minimal model whose 
results could be compared to data from a real system where the hypotheses 
are reasonably met [2H1 [25] ■ 

As stated previously, we seek to analytically solve the model introduced 
in [23] and test the consistence of the results obtained in [25]. Therefore we 
will be working on the set of differential equations that the authors derived 
for the dynamics of the fraction of monolingual speakers of each language. 
The derivation is as follows: 

We formulate the probability of shifting languages Pxy and Pyx and 
the probability of arriving to (departing from) the bilingual group from each 
monolingual group Pxb, Pyb (Pbx, Pby) based on the hypothesis listed 
above: 

Pxb = ck(l-s)(l-x) a , 
P YB = cks(l-y) a , 
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Pbx — Pyx 
Pby = Pxy 



c(l-k)s(l-y) a , 
c(l-k)(l-s)(l-x) a . 



(1) 



c > is a normalization constant [25] . Equations [T] are used to reckon the 
rates at which the population of the three groups grow or decline: 



A third differential equation exists that tracks the evolution of the proportion 
b of bilinguals, but this will not be needed in the following thanks to the 
normalization of the population x + y + b — 1. In [25] these equations are 
written such that the contributions of all terms from equations [1] can be 
explicitly read, but here we prefer a more compact notation so that we can 
research the field F = (F x , F y ): 



F x {x,y) = c[{l-x){l-k)s{l-y) a -x{l-s){l-x) a ], 
F y (x,y) = ci(l-y)(l-k)(l-s)(l-xy-ys(l- y y}. (3) 



For the present study, the parameters are restricted to k 6 (0, 1), s e (0, 1), 
and a > 0. Sometimes a will be further restricted. If this is the case, it will 



As said before, the parameter c normalizes the dynamics and is irrele- 
vant for equilibrium points and stability issues, therefore it was not paid 
much attention in previous literature and neither will it be paid attention 
now. The parameter a generalizes the monotonously increasing dependence 
of the probability of transition between languages as outlined in the fourth 
hypothesis. This parameter has been found to be larger than 1 and relatively 
constant among cultures (a ~ 1.31) in experimental accounts of the problem 
of language dynamics [T31 E2J E5] , but we intend to address the behavior of 
the system for a > 0, which is a more interesting generalization. The other 
parameters that concern us are k and s. 

3 Resolution of the model 

3.1 The model yields realistic trajectories 

We note that all the possible distributions of speakers among the different 
groups can be represented in the x—y space. There, the condition x+y+b = 1 



dx 

~dt 
dy 

dt 



F y {x,y). 



(2) 



be noted. 
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defines a triangular set A = {(x,y), x > 0, y > 0, x + y < 1} upon which 
the fields F x (x,y) and F y (x,y) are acting. For the sake of basic consistency 
of the model, its solutions must be feasible; meaning that a negative number 
of individuals in any group should be forbidden: the dynamics must happen 
inside A for realistic systems. 

Lemma 3.1 Assume the parameter a > 0. The set A = {(x, y), x > 0, y > 
0, x + y < 1} is positive invariant. 

Proof Let us see that the field defined in equations [3] is directed inwards in 
the boundaries of A. 

1. If x = and y G [0, 1]: 

F x (0,y) = c(l-k)s(l-y) a >0, 

F y {0,y) = ci(l-y)(l-k)(l-s)-ys(l- y y}. (4) 

The first inequality implies that the field flows inwards A since F x (0, y) > 
for < y < 1. 

2. If y = and x G [0, 1]: 

F x (x, 0) = c[(l-x)(l-k)s-x(l-s)(l-x) a ], 
F y (x,0) = c(l -k)(l -s)(l-x) a > 0. (5) 

For y = the field flows inwards A since F y (x, 0) > for any < x < 1. 

3. If y = 1 — x: 

F x (x,l-x) = c[(l- x)(l- k)sx a - x(l- s)(l- a;) a ], 
F B (a:,l-x) = c[x(l-A;)(l-s)(l-x) a -(l-x)s2; a ]. (6) 

In this case the field flows inwards A if and only if —F x (x, 1 — x) — 
F y (x, 1 — x) > since (—1, —1) is a normal vector to the straight line 
y = 1 — x pointing towards the interior of A. That condition is trivially 
satisfied and we conclude again that the field flows inwards into A 
through the segment y = 1 — x, x G (0, 1). 

This lemma means that any real distribution of speakers between the 
available groups that would evolve according to the proposed equations would 
remain feasible all the time. 
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3.2 Number of fixed points of the dynamics in A 

It is possible to find upper and lower limits to the number of equilibrium 
points that the system displays for different k, s, and a. We find the equilib- 
rium point of the system wherever the nullclines (curves defined by dx/dt — 
and dy/dt = 0) intersect each other. In A, the equilibrium points that can be 
detected by a simple inspection of the system are P x = (1, 0) and P y = (0, 1). 
We will term them trivial fixed points. There is another trivial fixed point for 
the dynamics, but it lays outside A: the point (1, 1). Depending on different 
values of the mentioned parameters we shall find more equilibrium points 
inside A. Following the notation introduced in this paragraph, we appreciate 
that the curves x = 1 and y = 1 are branches of the nullclines of the system. 
We name them the trivial branches. Our analysis will deal mainly with the 
non trivial branches. 

With a preliminary analysis of the nullclines of the field (F x = 0, F y = 0) 
it is possible to narrow down the number of fixed points in the interior of 
A to a maximum of 3: Equilibrium points of the dynamics are found in the 
intersections of the nullclines. Equating both components of the field to zero 
we get: 

(l-y) a = l ~ 8 1 x(l-x) a -\ 
s 1 — k 

Multiplying these equations and isolating y: 

y= n l , ■ (8) 

1 (l-k) 2 {l-x) 

We must restrict ourselves to x ^ 1 ^ y now to avoid divergences here and 
in following equations, but this is enought to continue with our discussion. 
Substituting equation [8] into the second expression of equations [3 

= 5(1 -*)(!-«) 



1 S)X \ 1 + ( 1 _ fe )2( 1 _ a .) 



— y \x+(i-kni-x)y=-^-(i-k). (9) 

x J 1 — s 



If (x*,y*) is an equilibrium point, x* must obey equation [9] and the corre- 
sponding y* is obtained from equation EJ 
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From the left-hand side of equation 



i \ a— 1 
1 — x 



g(x) = { — ) (x + (l-ky(l-x)) a . (10) 



If a > 1 it is true that: 



g(x) > 0, VxG (0,1); 
lim g(x) = +oo; 

x^0+ 

lim g(x) = 0. (11) 

Furthermore, g(x) has got a relative minimum and a relative maximum re- 
spectively at: 



Because all of this, the equation g(x) = ^£-(1 — k) can only have one, two, 
or three solutions in x G (0, 1) for fixed k and s, restricting thus the number 
of equilibrium points of the whole system. 

For a — 1, g(x) reduces to a straight line that might or might not fulfill 
g(x) = — k) within the range of interest x € (0,1). Because g(x) is 

a straight line, this equality can be obeyed for just one value of x at most. 
Thus for a = 1 there is at most one more fixed point within A, but it must 
not necessarily exist. 

Finally, for a < 1 the limits found in equation [TT1 swap: 

lim g(x) = 0, 

lim g{x) = +oo; (13) 

x— s>l~ 

and g{x) > in the whole range x G (0, 1). It is also monotonically increasing 
within this range and thus must always match j^-(l — k) in exactly one point 
x* internal to A. So for a < 1 there is always one fixed point besides the 
trivial ones. 
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3.3 Stability of the equilibrium points P x = (0, 1) and 

P y = (1,0) for a > 1 



We assess the stability of the system by evaluating the matrix at the existing 
fixed points, diagonalizing it, and considering the sign of the eigenvalues. For 
non trivial equilibrium points it becomes complicated to exactly locate them 
on the x — y plane, left asside its analysis through the Hessian matrix; but 
for P x and P y and restricting ourselves to a > 1 we can evaluate the Hessian 
matrix explicitly and it happens to be diagonal already: 



Furthermore, the eigenvalues are negative meaning that P x and P y are asymp- 
totically stable independently of the values of k and s for a > 1. 

Because P x and P y are always stable for a > 1 and the field is such that all 
trajectories enter A, if there is only one more equilibrium point x* interior to 
A it must be a saddle point and lie exactly at the frontier between the basins 
of attraction of (0, 1) and (1,0). If x* were unstable yet not a saddle point, 
either there would exist two more fixed points where the boundaries between 
basins cross the frontier of A, or there would exist trajectories leaving A; 
and neither of these is the case. If x* were stable it would have a basin of 
attraction for itself and new fixed points would need to exist in the separation 
between different basins. 

3.4 Studying the field in different regions of A for a > 1 

Now we will get more insights about the dynamics by further characterizing 
the field F at the boundary of A and in its interior. For this analysis we 
must assume a > 1, otherwise some of the functions that we will be making 
use of will be ill-defined. 

3.4.1 Studying F at the boundary of A 

Recalling F y (0,y) from equation HI we introduce: 




(14) 




(15) 
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and we note that it is continuous on the interval [0, 1), strictly decreasing on 
(0, 1/a) and strictly increasing on (I /a, 1). Since: 

G„(0) = G y (l) = (l-k)(l-s)>0, (16) 

we can find out if this function ever changes its sign by evaluating it at its 
minimum: G y (l/a). We get either: 

which would imply that F y (0,y) > Wy G [0, 1); or: 

(a - ir "- > (18) 



a" 

which would imply that there would exist yi, y 2 G (0, 1) such that F y (0, y) > 
if y G [0,yi) U (y 2 , 1] and F y (0,y) < for all y G (yi,y 2 ). In this case F y is 
zero at (0,yi) and (O,^)- Let us note that y% = y 2 if the equality holds on 
equation [181 

Likewise, recalling F x (x, 0) from equation [51 we define: 

G x (x) = (1 - k)s-x(l -s)(l -s) - 1 , (19) 

which is continuous on a; G [0, 1), strictly decreasing on (0, 1/a) and strictly 
increasing on (1/a, 1). Also: 

G X (0) = G X (1) = (l-k)s>0, (20) 

thus we find either: 

(a-l)"- 1 s 



< z (1-*), (21) 

1 — 6' 



which would imply F x (x, 0) > Wx G [0, 1); or: 

(a-1)"- 1 s 



> ; (1-*), (22) 

u~ 1 — S 

which would imply that there would exist xi,x 2 G (0, 1) such that F x (x, 0) > 
if x G [0, xi) U (x 2 , 1] and -^(x, 0) < for all x G (xx, x 2 ). In this case F x 
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is zero at (xi,0) and (x 2 ,0). Once again: x\ = x 2 if the equality holds on 
equation 

Let us note that if the strict inequalities [T7] and [22] are simultaneously 
true then s < 1/2 and if the strict inequalities [T8l and I2T1 are simultaneously 
true then s > 1/2. 

On the diagonal x + y — 1 the field takes the form written in equation 
[6] and the signs of F x (x, 1 — x) and 1 — x) can be studied analyzing 

(— -)° and (yz^) a respectively: 

Since the function (^-^) a 1 is strictly decreasing on (0, 1) and: 



fixed k and s. Therefore F x (x, 1 — x) < if x G (0, 2^) and i^fx, 1 — x) > 
if x G 1). 



Then F y (a;, 1 — x) > if x G (0, ^) and ^(x, 1 — x) < if x G (z y , 1). 

It can be trivially shown that z y < z x . Also it is true that \F x (x, 1 — x)| > 
Fj^x, 1 — x) if x G (0, z y ) and F x .(x, 1 — x) < |i^(x, 1 — x)| if x G (z x , 1). 

3.4.2 Further study of F x 

If x 7^ and x 7^ 1, then the points (x,y) which nullify the first component 
of the field F are those that obey: 




(23) 



then there exists only one z x G (0, 1) such that ( j = (1 — k)j^ for 





y 




(24) 



The function h x : [0, 1] — > M is defined as: 




(25) 
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on x G (0, 1] and h x (0) = 1. It is strictly decreasing on (0, 1/a) and increasing 
on (1/a, 1), it has got a minimum at 1/a, and h x (0) = 1 = h x (l). 

If y < h x (x) then F x (x,y) > 0. If y = h x (x) then F x (x,y) = 0. If 
y > h x (x) then F x (x,y) < 0. 

If the parameters of the system are such that inequality [21] holds, then 
h x (l/a) > and the plot of h x (x) intersects the boundary of A at (0, 1) and 
(z x , 1 - z x ) (figs. [TJi-b). 

If the strict inequality 1221 holds true then h(l/a) < and the plot of h x (x) 
intersects the boundary of A at (0, 1), (xi,0), (x2,0), and (z x , 1 — z x ) (figs. 
Hb-e). 

Additionally, since F y (x 2 , h x (x2)) > and F y (z x , h x (z x )) < 0, because F y 
and h x are continuous, it is warranted the existence of p G (x2, z x ) such that 
F(p, h x {p)) = (0, 0) -i.e. a fixed point of F. 

3.4.3 Further study of F y 

If y and y 7^ 1, then the points (x, y) which nullify the second component 
of the field F are those which obey: 

We define h y (y) similarly as we defined h x (x). This function is strictly de- 
creasing on (0, 1/a) and increasing on (1/a, 1), and ^(0) = 1 = h y (l). 

If inequality [T71 holds true then h y (l/a) > and the curve {(h y (y) , y) , y G 
[0, 1]} intersects the boundary of A at (1, 0) and (z y , 1 — z y ) (fig. □>). 

If the strict inequality [T8l holds true then h y {l/a) < and the curve 
{(h y (y),y),y G [0,1]} intersects the boundary of A at (1,0), (0,?/i), (0,y 2 ), 
and (z y , 1 — z y ) (figs. [T]a, anddb-e). 

Additionally, since F x (h y (y2), y 2 ) > and F x (h y {l — z y ), 1 — %) < 0, be- 
cause i 7 ^ and a y are continuous, it is warranted the existence of g G (1/2, ^~ z y) 
such that F(h y (q), q) = (0, 0) -i.e. a fixed point of F. 

The evolution of both h x (x) and h y (y) as a function of the parameters a, 
s, and k is partially shown in figs. [l]and|2j and can be dynamically explored 
in [281. 
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3.5 The nature of the orbits help us assess the stability 
of non-trivial fixed points 

The nullclines are always landmarks of the dynamic system under research. 
Their obvious use is to locate the equilibrium points in their intersections, 
but more information can be extracted if we look at them carefully. In section 
13.41 we used them to find out how the field behaves in the boundaries of A as 
they mark the sets of points where the vertical and horizontal components of 
the field are nullified. This applies also in the interior of A: The trajectories 
of the system pass by with vertical tangent through the points of the curve 
F x (x,y) = 0, and with horizontal tangent through the points of the curve 
F y (x,y) = 0. But also, these curves divide A in regions within which the 
signs of F x and F y are well determined. Topological arguments regarding the 
action of F upon these different regions of A help us put some limits to the 
kind of orbits that the system can yield: we will see that periodic dynamics 
can be banned. These considerations also let us find out whether non-trivial 
points are stable or not for a > 1. 

In this range of a the system will always have at least one more equilibrium 
point in the interior of A. We have seen that this can be deduced either from 
the crossings of h x (x) and h y (y) with the boundary of A or from equation 
[9] attending to the shape of g(x). The analysis of g(x) let us further know 
that also two or at maximum three fixed points can exist inside A. These 
three, two, or one equilibrium points will show up depending on the values 
of the parameters s, k, and a. Many possibilities are illustrated in fig. [TJand 
in [28]. 

With this in mind, let us consider the following regions: 

R x = {(x, y) G A/F x (x, y) < 0, F y (x, y) > 0}, 

R 2 = {(x,y) E A/F x (x,y) >0,F y (x,y) <0}. (27) 

Equivalently: 

i?i = {(x,y) G A/x < h y (x),y > h x (x)}, 

R 2 = {(x,y) G A/x > h y (x),y < h x (x)}. (28) 

In a similar way we could introduce: 

B x = {(x,y)eA/F x (x,y)>0,F y (x,y)>0}, 

B 2 = {(x,y) e A/F x (x,y) <0,F y (x,y) <0}; (29) 
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but these will not be interesting for us right now. 

Focusing on Ri and R 2 , they have got one or two connected components 
depending on if inside A there exist one, two, or three equilibrium points. 
We shall write R\ = A\ U A 3 , R 2 = A 2 U v4 4 ; being A 3 or A4 empty if on 
the interior of A there are not three equilibrium points, and A\ and A 2 the 
regions whose boundaries contain respectively P y = (0,1) and P x = (1,0). 
An account of these regions for some values of the parameters can be seen in 

fig. El 

Taking into account the sign of the components of the field we can tell 
that regions A±, A 2 , A3, and A 4 are positive invariant. If (xo,yo) G A\ then 
its trajectory (x(t),y(t)) for t G I t [0, 00) lays in A x because in the boundary 
of Ai the field points inwards. This trajectory is thus contained in a compact 
for t G It- Also, since < and ^jp- > 0, it can be verified that x(t) 
is monotonously decreasing on I t and y(t) is monotonously increasing on I t . 
Consequently it exists the limit \im t ^ +O0 (x(t),y(t)) = P y . The set A\ is 
therefore contained in the basin of attraction of P y . Analogously, it can be 
shown that if (xo,yo) G A 2 then \im t ^ +oc (x(t),y(t)) = P x : its trajectory lays 
in the region of attraction of P x ; thus, the basin of attraction of P x contains 
A 2 . 

Also, if A 3 and A 4 are both non-empty and (x , yo) G A3UA1 its trajectory 
remains either inside A 3 or inside A4 and converges towards the equilibrium 
point at the intersection of the frontiers of these regions. Since A 3 and A^ are 
non-empty only when there are three equilibrium points inside A, this result 
means that one of these three points, whenever they exist, must be stable. 
In this case we can determine that the two remaining fixed points in the 
interior of A must be saddle points. We do so with an argument similar to 
the one we used to show that x* is a saddle point when only one equilibrium 
point exists inside A (section l3.3p . We further deduce that the case with two 
interior fixed points corresponds to a saddle-node bifurcation and that this 
situation is the frontier between those cases with one and three equilibrium 
points in the space of parameters k — s. 

Concerning the dynamics of the system, it is important the following 
lemma which ensures that there is no oscillatory behavior: 

Lemma 3.2 There are not any periodic orbits in the x-y plane. 

Proof Actually, because of the regularity of the field, applying the Poincare- 
Bendixson theorem it can be deduced that if there would exist any closed 
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orbit it must enclose a fixed point on its interior. Thus, the periodic orbit 
would necessarily enter and exit two of the regions A^. This cannot happen 
because all regions Ai are positive invariant. 

This same argument also implies that fixed points cannot be foci, because 
trajectories approaching them should cross many times the frontiers between 
regions, some of which are positive invariant and cannot be left. 

3.6 Tentative solutions for a < 1 

Splitting the problem in a > 1 on the one side and a < 1 on the other made 
its solution easier because several of the reasonings that work very well in 
the former case are built on functions that are ill-defined in the later. An 
example are the functions G x (x) and G y (y), but also the Hessian matrix 
in (1,0) and (0,1) present some problems for a < 1. Luckily enough, in 
section [3721 we proved that there is just one more equilibrium point (x*,y*) 
for a < 1 which always appears if a < 1 and that might not appear for a = 1 
depending on the parameters, so we do not need to investigate 3 prospective 
fixed points as for a > 1. Also it is still valid the demonstration that A is 
positive invariant made in section 13.11 

In fig. [2] the nullclines are represented for various values of a and fixed k 
ands. We can observe (x*,y*) in the intersections, and we can also observe 
how the nullclines suffer a deep transformation as values of a larger than 1 
are employed. We shall study now the cases a = 1 and a < 1. Because of 
the analytic results are not so satisfactory, we shall complement them using 
numerical simulations whenever it is useful. These results should be ques- 
tioned as long as a complete mathematical proof is not available. 

It is particularly illustrative the resolution of the stability of P x = (1,0) 
and P y = (0, 1) for a = 1. This can still be analytically done. The Hessian 
matrix in this case reads: 

np/n IN ( - C (!- S ) -C{l-k)s \ 

lU ' > ~ \ -c(l - k)(l - s) + cs J ' 

- (-t(i- ,S ,Ki ( - 1 s r ) I)- (30) 

The eigenvalues are Ap = — c(l — s) and Ap = — c(l — k)(l — s) + cs for 
P y and Ap^ = — c(l — k)s + c(l — s) and Ap^ = — cs for P x . We see that 
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one of the eigenvalues is always the sum of two terms with different sign and 
this compromises the stability of P x and P y . Indeed, their stability depends 
now on the parameters k and s. By equating the conflictive terms to zero we 
obtain two curves relating k and s: 

1- kp 

2- k P , ' 

y 

(31) 

2-k Px { } 

These curves tell us where does the stability of P x and P y change in the space 
of parameters k — s: P x is stable for s > sp x and P y is stable for s < sp y . 
There is a region of values sp y < s < sp x where neither P x nor P y are stable 
and, since there are not any trajectories leaving A, there must exist a point 
(x*,y*) interior to A that is stable. 

The curves Sp x {kp x ) and sp y {kp y ) are plotted in fig. Hh- There it is shown 
the stability of the different stable points for different values of a in the k — s 
space, but aided by computer simulations. We see that the numerical results 
match the analytical results for sp x {kp x ) and sp y (kp y ), and that these curves 
seem to evolve into the boundaries between different regimes as a takes val- 
ues larger than 1. 

For a < 1 it is not possible to work out the stability of P x nor P y in a 
rigorous way. The Hessian matrix has diverging terms in this case: it is not 
well defined. We know, though, that there is always a third fixed point inside 
A and computer simulations suggest that this interior equilibrium point is 
always stable for a < 1. This would be in agreement with results obtained 
for the simpler precursor model [13] where a similar stable point is found for 
a < 1 [201126]. 

4 Discussion and conclusions 

The analytical results that we reached here are in partial discordance with 
those reported from previous numerical works [25] . Although the number of 
fixed points reckoned now could agree with previous accounts, the stability 
of them has been misread because of reasons that will be obvious right now. 
Notwithstanding this, the conclusions from [25] remain largely the same, as 
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we will see. We first discuss the results for a > 1 to compare directly with 
previous works and then we make some remarks about the cases a < 1. 

In [25] the stability of the system was assessed through computer sim- 
ulations only, thus stable fixed points were partly identified. Non-stable or 
saddle points did not stand out in these simulations because specific tests 
were not run therefore. It was concluded that one, two, or three stable equi- 
librium points existed in A, including the trivial ones at its boundary: P x , 
P y , and P* = (x*,y*); the later being the only equilibrium point explicitly 
detected in the interior of A. This matches exactly the picture drawn from 
the current work. The discrepancies arise regarding the stability of P x and 
P y : The computational tools used in [23J yielded a result that strongly sug- 
gested that this stability depended upon k and s (dependence upon a was 
not addressed: it was taken a = 1.31 for historical reasons) and that P* 
was stable whenever it could be detected by simulations-as it could only be 
detected if it was an attractor of the dynamics in the discretized version of 
equations [2j 

A plot was elaborated in [23] that divided the k — s space of parameters 
in five regions that would correspond to five different stability/instability 
combinations of the equilibrium points. Namely: i) only P x or only P y is 
stable, P* is not detected; ii) both P x and P y are stable, P* is not detected; 
iii) the three possible points are detected and stable; iv) P* and either only 
P x or only P y are stable; and v) only P* is stable. In all five cases P x and P y 
were supposed to exist and to be instable whenever their basins of attraction 
were found empty by the computer simulations. An updated version of that 
plot is reproduced in fig. H]c with the same five regions colored with different 
shades of blue and green. 

The interpretation given in [23] was not right although it was consistent 
with the numerical outcome. For example: we now know that P x and P y 
are always stable for a > 1, disregarding the values of k and s. We also 
know that at least one equilibrium point exists always in the interior of A, 
which may not be stable and which may lay in the boundary between the 
basins of attraction of P x and P y . In fig. [I] they are shown the plots of the 
nullclines for the exact same parameters as those used in |25j and we readily 
see how for certain parameters some saddle points interior to A approach P x 
or P y leading to a reduction of their basins of attraction. We now know that 
these equilibrium points never collapse into an unstable point. A basin of 
attraction may become undetectable to numerical means, and thus P x and P y 
may be deemed instable; but we have now found out analytically that P x and 
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P y remain stable for any value of k and s if a > 1. We were able to reproduce 
this numerical effect for the parameters used in [25] and for many others, as 
it can be seen in fig. HJa-d: lighter shades of blue or green indicate sets of 
parameters for which the computer simulations led to a wrong interpretation 
of the stability of some of the fixed points. 

For a > 1, the updated, more correct picture is as follows: There are al- 
ways 2 stable trivial fixed points P x and P y and depending on the parameters 
k and s there are 1, 2, or 3 more fixed points in the interior of A. It can be 
shown that if three points exist, one of them (termed P*) must be stable. 
It can be argued that if there is only one fixed point inside A it must be a 
saddle point; and that if there are three, those equilibrium points different 
from P x , P y , and P* must be saddle points as well. It can be guessed that 
the situation with 2 equilibrium points inside A corresponds to a saddle-node 
bifurcation as we transit through the k — s space from a region with one to 
a region with three non-trivial fixed points. The existing analytical evidence 
and fig. HI that shows the results of refined numerical simulations, are con- 
sistent with this view; although those facts that were not analytically proven 
in section [3] must be taken with enough care. 

These equilibrium fixed points and their stability have got a direct in- 
terpretation for the phenomenon for which the model was developed in the 
first place: they determine whether two coexisting languages would remain 
alive together, or if one of them is going to take over and extinguish the 
other. Also, the nature of the fixed points reached by the dynamics deter- 
mines whether individual bilingualism can be a stable trait. Given a pair of 
languages X and Y that coexist with status sx = s and sy = 1 — s and 
interlinguistic similarity k in a society with a fixed value of a > 1, the model 
presents two well differentiated regions: 

1. Coexistence is unstable: one language ends up suppressing the other 
and the bilingual group. What language survives depends on the ini- 
tial distribution of speakers among monolinguals of each language and 
bilinguals. This case corresponds to only one equilibrium point- which 
turns out not to be stable-in the interior of A and is depicted in figs. 
[Jji-b. The regions in the k — s space where this happens are colored in 
blue in fig. [4}>d. 

2. Coexistence is possible depending on the parameters k, s, and a; 
and on the initial conditions of the dynamics. This case is the one with 
three stable equilibrium points P x , P y , and P*. The initial conditions 
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determine whether a language drives the other to extinction and makes 
bilingualism disappear, or if a steady state is reached (P*) in which 
groups of monolingual speakers of both languages survive along with 
a bilingual group. This is what happens in figs. H]c-e; and parameters 
k — s for which we find this situation are indicated in green in figs. 

Sfc-d. 

A third case regarding number of fixed points would exist at the boundary 
between these regions, but it does not seem to introduce any new behavior 
attending to the coexistence of languages. Coming back to the incomplete 
interpretation made of the results in [22], the two cases just outlined already 
include those configurations of parameters for which some attractors are so 
small that the extinction of a language or the coexistence of both of them 
is almost unavoidable without regard of the initial conditions; although we 
now know that this is never the case. 

We can see from the numerical simulations in fig. SI that regions of the 
k — s space where stable bilingualism is possible correspond to those with a 
more balanced status between languages. This balance is not so important 
for larger interlinguistic similarity: then a stable bilingual situation can be 
reached even for well distinct sx and sy, depending on the initial distribution 
of speakers. Further illustration of the role of k and s is made in fig. that 
shows qualitative bifurcation diagrams of the stable fixed points when varying 
these parameters with fixed a. 

The parameter a was found relatively constant among cultures as indi- 
cated before [33], and this justified why it was not payed that much attention. 
But now we have also studied how the possible outcomes change as a varies. 
First, considering only a > 1, we see (fig. H]) that a larger a means that 
the possibilities for stable bilingualism are reduced. This parameter was al- 
ready considered in an analysis of the more basic Abrams-Strogatz model 
[T7J [201 [26] , and it was cleverly termed volatility parameter: the lower a the 
more volatile a large group becomes and vice- versa. Thus, for larger a bigger 
groups are more persistent and it is smaller the set of parameters for which 
it can be reached a more diluted distribution of speakers (this would be: 
a solution with speakers belonging to the bilingual groups, or communities 
with monolingual groups of each language coexisting together). For lower 
a, larger monolingual groups are not so permanent and a steady solution is 
easier to reach in which all languages coexist. 

This volatility is a critical feature at a = 1: P x and P y become unstable 
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for some of the parameters k and s as it was said before. In terms of language 
dynamics this means that a monolingual group agglutinating all the speakers 
is no longer possible, whatever the initial conditions, if their statuses are close 
enough and depending on the similarity between languages. The region of 
the k — s space where this happens can be seen in green in fig. H^i. In 
such cases the only stable solution in the long term is the coexistence of the 
monolingual groups along with the bilingual one. But still at a = 1 one 
language might extinguish the other if its status is larger enough. There 
is a crucial difference between such an extinction and those happening for 
a > 1: before, both languages could survive depending on the initial number 
of speakers of each one; now the extinction does not depend on this initial 
condition if the parameters are those needed for a language to take over (blue 
regions in fig. Hk). 

The possibility of language extinction seems to change completely for 
a < 1: then the volatility is so high that the monolingual options are never 
stable and the survival of both languages within their monolingual groups 
and along a bilingual group of speakers is guaranteed for any values of the 
parameters k and s, and for any initial distribution of speakers. It was not 
possible to prove this very last result analytically beyond any doubt and it 
was obtained thanks to numerical simulations. This solution is consistent 
with similar outcomes for the seminal Abrams- Strogatz model [T7J [201 126] - 
which should be the limit case for k — > and b — > of the equations under 
research in this paper. 

4.1 Nature of the orbits and higher order contribu- 
tions 

A very important contribution of this paper is that brought in by lemma 
13.21 It is clear its mathematical meaning: because of the nature of the field 
(F x , F y ) around a fixed point it is not possible to find closed-i.e. periodic- 
orbits. The same lemma also implies that any solution must consist of an 
exponential decay towards a fixed point: that an oscillatory decay is not 
possible. The interpretation of this result is rather strong when it comes 
down to languages: the extinction or raise of languages must be a monotonous 
phenomenon according to the present equations. Tendencies that could be 
expected, e.g. alternation in the preponderance of a language in a region, 
should not be observed. If such result were derived from real data, research 
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should be focused on what is needed to complement the present model: "what 
would be the minimum elements that play a role in generating cyclic behavior 
in language competition?" , because those employed here would not suffice. 

This lemma has got also some implications even if we would consider 
higher order or stochastic extensions of the present model. The equations in- 
vestigated in this paper are nothing but a deterministic, mean field approach 
to a phenomenon that usually takes place on a stochastic environment. The 
next more realistic strategy to model language competition or coexistence 
departing from our current equations would be simulations of discrete agents 
that shift between the monolingual or the bilingual groups at random, be- 
ing the transition probabilities given by equations [TJ This is coherent if we 
suppose some free-will and variability to the speakers when deciding what 
language to use. We expect thus intrinsic stochasticity to be present and 
manifest throughout noise. The power spectrum of this noise can be investi- 
gated. The nature of the equilibrium points found for equations [2] establishes 
some important limitations to the kind of dynamics that can arise, as we will 
argue. 

Several techniques are available in the literature to incorporate uncer- 
tainty in a deterministic model, from agent-based simulations with stochastic 
interaction events [29j [30] to theoretical considerations of a more analytical 
nature [311 [32j [33]. According to |M] . we can determine that some inter- 
esting phenomena are ruled out from our model because of the lack of foci 
fixed points (again, recalling lemma I3~2"j) . Namely, it is not possible to find 
Stochastic Amplification of Fluctuations (SAF), a phenomenon that offers a 
possible explanation to emergent quasi-oscillations observed in fields as di- 
verse as ecology [33], epidemiology [35], or brain dynamics [Ml [36]. In SAF 
the spectrum of the noise would present a prominent peak corresponding 
to these quasi-oscillations. Opposed to this, in our study case the power 
spectrum of intrinsic noise must present a monotonous decay proportional 
to l/u 2 , and no outstanding peaks. We confirmed this result for many sets 
of parameters with agent-based simulations, as suggested above, finding no 
interesting features in the spectra, in agreement with the SAF theory. 

SAF, if present, would be self evident in agent-based simulations. It could 
also be noted in series of real data if reports of language usage over time with 
enough precision were available. SAF is usually associated with adaptive 
reacting forces such as prey-predator or activator-repressor dynamics. Thus, 
reports of emergent oscillations in language dynamics could warn us ofthe 
presence of such forces driving language competition and serve for further, 
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necessary refinement of the model. Numerical evidence shows that periodic 
solutions may appear if the status s of the languages were allowed to change 
over time, which is a rather realistic extension of the model. Also, if there 
would exist models of language coexistence that presented SAF in a natural 
way, the observation (or the not observation) of this phenomenon in real data 
could help us determine which one is closer to reality. 
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Figure 1: Non-trivial nullclines of the system for different values of 
k and s and fixed a = 1.31. Five pairs of parameters (k, s) were used 
to plot the non-trivial branches of the nullclines. The values were chosen 
to compare our results with those from [25J, thus: a k = 0.65, s = 0.80, 
b k = 0.20, s = 0.40, c k = 0.65, s = 0.50, d k = 0.75, s = 0.35, and 
e k = 0.75, s = 0.35. Stable fixed points are indicated in red and saddle 
points in green. Contrary to what was interpreted in [25], cases a and b on 
one side and c, d, e on the other side are equivalent to each other. The five 
differentiated cases are obvious, though, if we attend to the sizes of certain 
basins of attraction: some of them can hardly be detected with numerical 



systems. We can appreciate in the different figures how important features, 
like crossings with the boundary evolve for the different parameters. This 
can be better grasped in [28] . 
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Figure 2: Non-trivial nullclines of the system for different values of 

a. Fixed k = 0.65 and s = 0.50 is taken except in panel c. This figure 
illustrates how the non- trivial branches of the nullclines suffer a sudden 
change as the parameter a starting from a < 1 increases through a = 1 and 
to a > 1. a For a < 1 the nullclines must always cross inside A. b a = 1 is 
the only case when the nullclines are straight lines that go through (0, 1) and 
(1,0) respectively. As straight lines, depending on their slopes they might 
cross inside A, but this is not necessarily the case (c a = 1, k — 0.3, and 
s = 0.75). d For a > 1 the shape of the nullclines becomes more complicated 
and the possible crossings must be carefully addressed. We saw the many 
possibilities in fig. [TJ 
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Figure 3: Positive invariant regions Ri and R 2 on the x — y. We 

identify R\, where F x is negative and F y is always positive (green in the 
figure); and R 2 , where F x is positive and F y is negative (red regions), a 
a = 1.40, s = 0.50, k = 0.65. Three equilibrium points exist inside A and 
regions A3 and A4 are non-empty. Both R\ and R2 present two connected 
components (A\ and A3, and A 2 and A4 respectively) which are, each of 
them, positive invariant. This means that the dynamics do not exit any of 
these regions once they enter: they must tend to a stable fixed point in their 
boundary. Thus, we see how points inside A 1 are taken to P y and points 
inside A 2 are taken towards P x . The only possibility for regions A 3 and A 4 is 
that it exists another stable fixed points exactly in the joint between the two 
of them, b a = 1.31, k = 0.20, s = 0.40. In this case only one equilibrium 
point exists inside A. Regions A3 and A4 are empty, but the same as before 
applies to regions A\ and A 2 . points in their interior must be driven towards 
P y and P x respectively. 
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Interlinguistic Similarity (k) 

Figure 4: Regions of the k — s space with different number of equi- 
librium points, a a=l: It is possible to compute analytically the curves 
sp y (kp y ) (thick black line) and sp x (kp x ) (thick red), which determine the fron- 
tiers above (below) which P y (respectively P x ) are unstable. The green area 
in which both P y and P x are unstable need the existence of an equilibrium 
point inside A which is stable to attract the dynamics. The blue regions 
correspond to values of the parameters where either P x or P y are stable, 
but only one of them; meaning, in terms of competing languages, that one 
tongue must extinguish the other whatever the initial conditions. While the 
curves sp y {kp y ) and sp x (kp x ) were computed analytically, the colored regions 
were found out through computer simulations: a point of the k — s space 
would be painted in blue if either of the attractors P x or P y were found after 
evolving the system a time large enough, and it would be painted in green 
whenever the simulations did not converge towards P x nor P y in a similar 
umber of iterations. Numerical and analytical results agree, b a = 1.1, c 
a = 1.31, da = 1.5: In either of these cases blue regions indicate that only 
the attractors P x and/or P y have beeflOietected: both of them were detected 
in dark blue regions while only one of them was detected in lighter blue re- 
gions. Green regions indicate now that an attractor P* interior to A has been 
detected: the darkest area correspond with P* being detected along P x and 
P y and the lighter areas means that one or two of P x and P y have not been 
detected. All these attractors were found numerically. The different shades 
of blue and green aim at demonstrating how the stability of the fixed points 




Figure 5: Bifurcation diagrams. Qualitative bifurcation diagrams are 
shown for fixed a = 1.31 with: a varying k and fixed s = 0.5 and b varying 
s and arbitrary k close to 0.8. Beginning with low k, when varying this 
parameter while holding s fixed we depart from a situation in which P x and 
P y are stable and only one extra equilibrium point exists in the interior of 
A, a saddle point. Above a certain value of k a new stable fixed point P* 
comes into existence. For varying s and fixed k, if k is large enough-as it is 
the case in this example-they exist two values of s between which the third 
stable equilibrium P* point exists. The trajectory of P* through the x — y 
space is qualitative in both cases. 
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